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1.  SUMMARY 


The  effort  described  herein  seeks  to  explore  the  near-resonant  thermomechanics  of  energetic  and 
mock  energetic  particulate  composite  materials.  The  effort  specifically  focuses  on:  (i) 
characterizing  the  macroscale,  elastic  and  plastic  responses  of  these  materials  under  various 
mechanical  excitations  at  a  range  of  ambient  temperatures;  and  (ii)  developing  preliminary 
computational  modeling  tools,  which  can  be  used  to  predict  material  response  during  energetic 
material  formulation  and  munition  design.  Key  topics  described  herein  include:  sample 
preparation;  macroscale  thermomechanical  modeling  and  experimentation;  dissipative  modeling 
and  system  identification;  endochronic  plasticity  modeling;  and  crystal-binder  interface 
modeling. 

2.  INTRODUCTION 

The  effort  described  herein  seeks  to  explore  the  near-resonant  thermomechanics  of  energetic  and 
mock  energetic  particulate  composite  materials,  building  upon  both  the  prior  work  of  the  Pis 
(see,  for  example,  [1-6])  and  prior,  related  work  in  the  field  pertaining  to  the  periodic  excitation 
of  energetic  materials  (see,  for  example,  [7-8])  and  hot-spot  formation  (see,  for  example,  [9]). 
The  effort  specifically  focuses  on:  (i)  characterizing  the  macroscale,  elastic  and  plastic  responses 
of  these  materials  under  various  mechanical  excitations  at  a  range  of  ambient  temperatures;  and 
(ii)  developing  preliminary  computational  modeling  tools,  which  can  be  used  to  predict  material 
response  during  energetic  material  formulation  and  munition  design.  More  specifically,  Task 
Order  0001  spans  three  inter-related  research  tasks: 

•  Task  1  -  Macroscale  Structural  and  Mechanical  Modeling:  This  task  emphasizes  further 
development  of  macroscale,  distributed-parameter  and  lumped-mass,  combined  structural 
and  thermal  models  of  mock  energetic  and  energetic  materials  subjected  to  mechanical 
vibration.  The  models  are  designed  to  be  amenable  to  a  variety  of  particle/binder 
material  systems,  a  range  of  input  excitations,  and  various  ambient  temperatures. 

•  Task  2  -  Modeling  of  Damping  and  Dissipation:  This  task  emphasizes  the  development 
of  refined,  macroscale  dissipation/damping  models  for  mock  energetic  and  energetic 
materials  subjected  to  mechanical  vibration. 

•  Task  3  -  Computational  Mechanics  Tool  Development:  This  task  emphasizes  the 
development  of  constitutive  models  for  mock  energetic  and  energetic  materials  that 
account  for  nonlinear  stress-strain  dependencies  due  to  damage,  rate-dependent  and  rate- 
independent  dissipative  properties,  as  well  as  temperature-dependent  properties,  and  their 
distillation  into  computational  tools  suitable  for  use  in  energetic  material  formulation  or 
munition  design. 

This  final  report  details  the  various  advancements  made  by  the  principal  investigators  and  their 
research  assistants  over  the  noted  period  of  performance.  The  reader  should  note  that  all  of  the 
research  activities  described  herein  are  continuing,  in  a  big-picture  sense ,  under  Task  Order 
0002,  and  thus  the  methodologies,  results,  and  perspectives  detailed  herein  represent  only  a 
snapshot  of  ongoing  research. 
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The  subsequent  sections  are  organized  in  a  traditional  sense  with  methodologies,  key  results,  and 
discussion  provided  in  turn.  Discussions  of  sample  preparation;  macroscale  thermomechanical 
modeling  and  experimentation;  dissipative  modeling  and  system  identification;  endochronic 
plasticity  modeling;  and  crystal-binder  interface  modeling  are  provided  in  Sections  3  and  4. 
Please  note  that  detailed  derivations  related  to  the  dissipative  modeling  and  system  identification, 
as  well  as  the  endochronic  plasticity  modeling,  are  presented  in  the  Appendix  to  enhance  the 
flow  the  report. 

3.  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 
3.1  Sample  Preparation 

Mock  energetic  materials  in  a  variety  of  geometries  have  been  fabricated  with  the  goal  of 
studying  the  effect  of  formulation  on  bulk  heating  due  to  mechanical  vibration.  The  formulation 
variations  are  based  upon  the  PBXN-109  formulation  given  by  Hamshere  et  al.  [10],  with  sugar 
replacing  the  RDX  unless  otherwise  noted.  A  summary  of  this  base  formulation  is  given  in  Table 
1. 


Table  1.  PBXN-109  Formula! 

tion 

Constituent 

Weight  Percent 

RDX 

64.00% 

Aluminum 

20.00% 

R45-HT 

7.346% 

(Hydroxyl-terminated  Polybutadiene  Resin) 
Dioctyl  Adipate  (DOA) 

7.346% 

Antioxidant  2246 

0.100% 

Dantocol  DHE 

0.260% 

Triphenylbismuth 

0.020% 

Isophorone  Diisocyanate  (IPDI) 

0.950% 

The  base  formulation  given  above  was  varied  based  on  two  parameters,  as  defined  below.  Each 
parameter  is  being  studied  at  three  levels,  yielding  a  sample  matrix  of  seven  possible  mock 
energetic  formulations. 

1.  Solids  Loading  (a  metric  to  indicate  the  amount  of  the  mixture  comprised  of  "solid" 
powder):  Levels  -  0%,  50%,  and  85%. 

2.  Additive  Content  (a  metric  to  indicate  the  amount  of  powder  in  the  mixture  is  comprised 
of  aluminum):  Levels  -  0%,  15%,  and  30%. 

The  formulations  are  cast  into  three  mold  geometries.  Both  small  and  large  cylinders  are 
prepared  for  uniaxial  compression  tests  in  order  to  determine  material  properties  and  model 
behavior.  Plate  geometries  are  prepared  for  harmonic  base  excitation  tests  to  investigate  thermal 
and  mechanical  response  at  the  bulk  scale.  The  mold  assemblies  used  are  shown  in  Ligures  1-3. 
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Figure  1.  Mold  assembly  for  the  2.54  cm  diameter,  2.54  cm  height  cylindrical  samples 


Figure  2.  Mold  assembly  for  the  7.62  cm  diameter,  7.62  cm  height  cylindrical  samples 
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Figure  3.  Mold  assembly  for  the  17.78  cm  x  25.4  cm  x  1.27  cm  rectangular  plate  samples 


Two  samples  of  each  geometry  are  cast,  yielding  a  sample  set  of  6  samples  per  formulation.  The 
samples  successfully  fabricated  at  this  time  are  summarized  in  Table  2.  The  50%  solids  loading 
samples  are  slated  to  be  completed  in  the  coming  months.  An  issue  of  solids  settling  out  of  the 
binder  during  curing  must  be  addressed  at  the  lower  solids  loading  before  fabrication.  Current 
research  seeks  to  "pre-cure"  the  HTPB  binder  to  increase  its  effective  viscosity,  allowing  it  to 
hold  the  solids  in  suspension  long  enough  for  solidification  to  occur. 

Table  2.  Sample  Preparation  Record  as  of  September  27,  2016 


Formulation 

Samples  Cured 

0%  Solids  Loading 
0%  Additive  Content 

(2)  2.54  cm  diameter/height  cylinders 
(2)  7.62  cm  diameter/height  cylinders 
(2)  17.78  cm  x  25.4  cm  x  1.27  cm  plates 

50%  Solids  Loading 
0%  Additive  Content 

(1)  2.54  cm  diameter/height  cylinder 
(1)  7.62  cm  diameter/height  cylinder 
(1)  17.78  cm  x  25.4  cm  x  1.27  cm  plate 

85%  Solids  Loading 
0%  Additive  Content 

(2)  2.54  cm  diameter/height  cylinders 
(2)  7.62  cm  diameter/height  cylinders 
(2)  17.78  cm  x  25.4  cm  x  1.27  cm  plates 

85%  Solids  Loading 
15%  Additive  Content 

(2)  2.54  cm  diameter/height  cylinders 
(2)  7.62  cm  diameter/height  cylinders 
(2)  17.78  cm  x  25.4  cm  x  1.27  cm  plates 

85%  Solids  Loading 
30%  Additive  Content 

(2)  2.54  cm  diameter/height  cylinders 
(2)  7.62  cm  diameter/height  cylinders 
(2)  17.78  cm  x  25.4  cm  x  1.27  cm  plates 
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In  regards  to  the  fabrication  process,  binder  constituents  are  mixed  manually  until  coarsely 
combined.  A  Resodyn  acoustic  mixer  is  then  employed  to  develop  a  nominally  homogeneous 
mixture  (samples  are  mixed  for  5  min  at  a  constant  intensity  of  80%).  The  solids  are  then  added 
into  the  mixture,  and  the  Resodyn  process  is  repeated  once  again.  The  formulation  is  then 
transferred  into  the  molds  described  above  and  cured  in  an  oven  at  60  °C  for  7  days. 

3.2  Macroscale  Thermomechanical  Testing 

Numerous  efforts  have  been  made  to  increase  the  understanding  of  particulate  composite  mock 
energetic  plates  under  mechanical  vibration.  As  part  of  this  effort,  experiments  were  conducted 
on  previously-fabricated  plate  samples  in  order  to  investigate  the  effect  of  thermal  boundary 
conditions  on  the  thermal  response  of  representative  samples.  The  investigation  into  the  effect  of 
formulation  variation  on  the  thermal  and  mechanical  response  has  also  begun,  with  the  results 
gathered  to  date  being  described  in  the  following  sections. 

As  alluded  to  above,  harmonic  base  excitation  experiments  were  performed  in  order  to 
investigate  the  effect  of  thermal  boundary  conditions  on  the  vibration-based  heating  of 
particulate  composite  plates.  The  plate  samples  tested  consisted  of  HTPB  with  varying  ratios  of 
embedded  ammonium  chloride  crystal,  which  were  chosen  as  a  rough  mechanical  mock  for 
ammonium  perchlorate  (AP).  Details  on  these  samples  can  be  found  in  References  [3]  and  [4]. 
These  samples  were  tested  under  both  ambient  and  insulated  environmental  conditions. 

Ambient  experiments  were  conducted  using  three  major  pieces  of  equipment.  The  plate  samples 
were  mounted  in  a  pre-existing  fixture  to  approximate  a  clamped-free-clamped-free  (CFCF) 
boundary  configuration  and  placed  upon  a  TIRA  59335/LS  AIT-440  electrodynamic  shaker.  A 
Polytec  PSV-400  scanning  laser  Doppler  vibrometer  was  then  employed  to  record  frequency 
responses  and  operational  deflection  shapes.  Simultaneously,  a  FLIR  A325  infrared  camera  was 
used  to  capture  the  temperature  distribution  of  the  top  surface  of  the  plate  using  infrared 
thermography.  The  experimental  setup  is  depicted  in  Figure  4.  The  results  obtained  from  this 
and  related  testing  can  be  found  in  Section  4.1. 
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Figure  4.  Labeled  schematic  of  the  experimental  setup  used  under  ambient  environmental 
conditions 


3.3  Dissipative  Modeling  and  Material  Property  Parameter  Estimation 

The  PBXN-109  samples  described  above  are  being  tested  under  this  portion  of  the  described 
research  effort  with  an  aim  towards  dissipative  modeling  and  material  property  parameter 
estimation.  Both  random  and  harmonic  base  excitation  tests  have  been  conducted  on  the  25.4  and 
76.2  mm  (1  and  3  in.)  tall/diameter  samples  with  excitations  over  the  frequency  range  from  10  to 
1000  Hz.  Base  excitation  is  provided  by  the  aforementioned  TIRA  electrodynamic  shaker  and 
the  samples  are  loaded  with  a  mass  and  constrained  to  move  uni-axially.  The  geometry, 
composition,  and  cure-date  of  the  samples,  as  well  as  a  testing  progress  summary  are  included  in 
Table  A1  in  the  Appendix. 

It  should  be  noted  that  a  day  of  harmonic  base  excitation  testing  consists  of  3  sweeps  through 
resonance  30  min  apart  at  1  g,  2  g,  3  g,  or  5  g.  Currently  we  are  conducting  two  excitation  level 
tests  on  the  same  day  and  the  tests  are  ordered  as  low,  low,  high,  high,  low,  high  (where  low 
denotes  the  lower  excitation  level  and  high  the  higher).  This  will  allow  us  to  examine  same  day 
repeatability  and  observe  if  the  order  of  the  base  excitation  level  affects  the  mass-material  system 
response.  A  day  of  random  base  excitation  consists  of  3  tests  with  excitation  over  the  range  from 
10  Hz  to  1000  Hz.  Testing  levels  are  specified  by  the  power  spectral  density  level  in  the 
frequency  range  of  the  excitation.  Levels  being  used  are  0.015  g2/Hz,  0.025g2/Hz,  or  0.050 
g2/Hz.  The  duration  of  each  test  is  100  s.  Note  that  in  the  course  of  testing,  some  excitation 
levels  had  to  be  adjusted  to  avoid  samples  fracturing.  The  mass  loading  of  the  samples  was  also 
adjusted  so  that  the  static  strain  levels  were  similar  for  both  sizes  of  samples. 
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Note  that  details  of  the  models  and  the  iterative  system  identification  process  that  is  used  to 
extract  the  material  parameters  from  the  experiments  described  in  this  section  are  provided  in  the 
Appendix. 

3.4  Endochronic  Plasticity  Modeling 

Energetic  composite  materials  may  exhibit  significant  nonlinear  stress-strain  response  under 
load,  without  a  distinctive  yield  surface.  As  a  result,  linear  material  assumptions  and  classical 
plasticity  theories  may  prove  to  be  insufficient  in  predicting  their  mechanical  response, 
necessitating  the  use  of  a  plasticity  theory  that  does  not  include  the  notion  of  a  yield  surface.  The 
endochronic  theory  of  plasticity  developed  by  Valanis  [11]  is  a  hereditary  plasticity  theory  that 
falls  into  this  category.  The  theory  is  hereditary  in  the  sense  that  that  the  state  of  stress  in  the 
neighborhood  of  a  plastic  material  point  depends  on  the  strain  or  deformation  history  of  that 
point.  The  strain  history  is  represented  by  a  variable  called  intrinsic  time,  and  the  stress  evolution 
is  given  by  the  convolution  integral  between  the  strain  tensor  and  the  memory  kernel,  which  is  a 
scalar  function  of  intrinsic  time.  For  a  plastically  incompressible,  rate-independent,  and  initially 
isotropic  material,  the  endochronic  constitutive  equation  is  given  by 


where  z  is  called  the  intrinsic  time  scale,  p(z)  is  the  kernel  function  and  s  is  the  deviatoric  stress 
given  by 


1 

s  —  a  —  -tr(a)I  (2) 

and  ep  is  the  plastic  strain.  The  total  strain  e  is  given  by 

e  =  ee  +  ep  (3) 

where  ee  the  elastic  strain,  related  to  the  stress  a  by  the  elastic  stress-strain  relationship 

a  —  X  tr(ee)  +  2 \iee  (4) 

where  X  and  p  are  the  Lame’s  constants. 

The  intrinsic  time  scale  z  is  related  to  another  variable  (  and  a  function  f(Q  by  the  relation 


dz  — 


m 


(5) 


where  (  is  called  the  intrinsic  time  measure  and  is  a  representative  of  the  plastic  strain  history  in 
the  form 


d (  =  yJdep\  dep 
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(6) 


and  /(£)  is  called  the  hardening  function  and  represents  the  hardening  or  softening  behavior  of 
the  material.  Thus,  the  intrinsic  time  scale  z  incorporates  both  strain  history  and  material 
hardening  or  softening.  Researchers  have  successfully  used  this  theory  to  model  the  inelastic 
behavior  of  many  materials  like  metals  [12,13],  concrete  [14],  soils  [15],  metal  matrix 
composites  [16],  filled  rubber  [17],  Asphalt  [18],  etc.,  which  confirms  the  validity  of  the  model 
and  its  applicability  to  our  problem.  Therefore,  endochronic  theory  has  been  chosen  as  the  basis 
of  our  predictive  computational  tool  for  energetic  composite  materials. 

In  order  to  formulate  a  computationally-efficient  and  robust  tool,  the  endochronic  theory  has 
been  implemented  in  the  form  of  a  numerical  scheme  reported  by  Hsu  et  al.  [19],  which  consists 
of  numerical  integration  of  the  endochronic  constitutive  equations  in  3D.  The  scheme  has  been 
modified  to  include  rate-dependent  (viscoplastic)  effects.  In  addition,  an  efficient  staggered 
algorithm  has  been  developed  for  simultaneous  calibration  of  the  hardening  and  kernel  functions. 
The  complete  numerical  scheme  along  with  calibration  procedure  has  been  described  in  the 
Appendix. 

3.5  Crystal-Binder  Interface  Modeling 

There  are  well-established  classical  experimental  and  analytical  studies,  which  explore  the 
bonding  between  polymers  and  rigid  particles.  However,  most  numerical  studies  of  fracture  in 
composite  materials  use  particle/matrix  interface  cohesive  laws  and  parameters  extracted  from 
the  cohesive  energy  at  the  macroscale  level.  While  such  an  approach  enables  good  results  when 
the  mechanical  response  is  compared  to  macroscale  experiments,  in  high-energy  composites 
where  the  localization  of  energy  dissipation  by  delamination  and  friction  may  lead  to  the 
initiation  of  hot  spots,  the  local  particle/matrix  cohesive  law  is  of  extreme  importance.  In  this 
work,  a  systematic  numerical  study  is  performed  to  establish  a  cohesive  law  of  the  composite 
components  and  interfaces  and  to  obtain  the  fracture  energy. 

Cohesive  zone  models  have  been  very  successful  in  simulating  fracture  in  many  types  of 
composite  materials  [20],  especially  when  cracks  nucleate  and  propagate  along  the  interfaces 
[21,22].  In  contrast,  phase  field  descriptions  of  damage  have  an  advantage  in  problems  where  the 
crack  path  is  unknown  a  priori  [23,24].  PBXs  are  a  special  class  of  composites  in  which  matrix, 
particles,  and  interfaces  are  prone  to  fracture,  and  therefore  a  phase  field  approach  is  more 
effective.  The  purpose  of  this  work  was  to  determine  the  cohesive  law  for  particle,  matrix  and 
interfaces  using  a  phase  field  damage  mechanics  (PFDM)  approach.  Furthermore,  the  results 
were  used  to  calibrate  the  material  parameters  needed  for  future  more  complex  simulations. 

PFDM  simulations  were  performed  to  predict  debonding  of  the  matrix-particle  interface  in  a 
spherical  glass  particle  embedded  in  Sylgard.  The  simulations  are  performed  using  a  3D  parallel 
software  (Note:  The  code  is  available  at  http://c-primed.github.io/fvm/),  MEMOSA,  which  is 
based  on  the  cell-centered,  finite- volume  method  [25].  The  number  of  cells  used  in  the 
simulations  is  123,048  and  they  have  an  average  size  20  pm.  The  particle  matrix  interface  is 
defined  as  a  region  with  twenty  elements  with  length  2  pm  in  the  radial  direction.  Details  of  the 
model  and  algorithm  can  be  found  in  Xie  et  al.  [26]. 
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4.  RESULTS  AND  DISCUSSION 


4.1  Macroscale  Thermomechanical  Testing  Results 

Although  the  HTPB/ammonium  chloride  plate  samples  had  been  previously  tested  under  ambient 
conditions  in  prior  work,  the  tests  were  repeated  in  order  to  gather  baseline  data  and  quantify  the 
impact  of  sample  aging.  The  vibrometer-recorded  HI  frequency  response  estimators  for  a 
representative  plate  sample  in  ambient  conditions  in  response  to  three  levels  of  excitation  are 
presented  in  Figures  5-7  for  the  50%,  65%,  and  75%  solids  loading  plates,  respectively. 


Frequency  (Hz) 

Figure  5.  Experimental  HI  mechanical  frequency  response  estimator  for  a  50%  solids 
loading  plate  obtained  at  three  levels  of  excitation.  The  blue,  green,  and  red 
curves  depict  responses  at  1, 1.86,  and  2.44  g  RMS,  respectively.  Solid  lines 
represent  data  from  the  geometric  center,  and  dashed  lines  represent  data 
from  the  offset  point. 
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Figure  6.  Experimental  HI  mechanical  frequency  response  estimator  for  a  65  %  solids 
loading  plate  obtained  at  three  levels  of  excitation.  The  blue,  green,  and  red 
curves  depict  responses  at  1, 1.86,  and  2.44  g  RMS,  respectively.  Solid  lines 
represent  data  from  the  geometric  center,  and  dashed  lines  represent  data 
from  the  offset  point. 


Frequency  (Hz) 

Figure  7.  Experimental  HI  mechanical  frequency  response  estimator  for  a  75%  solids 
loading  plate  obtained  at  three  levels  of  excitation.  The  blue,  green,  and  red 
curves  depict  responses  at  1, 1.86,  and  2.44  g  RMS,  respectively.  Solid  lines 
represent  data  from  the  geometric  center,  and  dashed  lines  represent  data 
from  the  offset  point. 
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Concurrent  with  mechanical  excitation  and  characterization,  the  thermal  response  for  each  of  the 
plate  samples  was  recorded  using  the  aforementioned  FLIR  infrared  camera.  These  experiments 
were  performed  over  a  60  min  time  interval  under  a  2  g  sinusoidal  excitation  near  the  first 
resonant  frequency  of  each  respective  plate  (where  the  highest  heating  is  typically  expected).  The 
measured  plate  surface  temperatures  are  plotted  in  Figures  8  and  9.  The  legend  indicates  the 
solids  loading,  followed  by  the  sample  number.  The  colored  envelope  surrounding  each  curve 
indicates  one  standard  deviation  for  each  trial. 
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Figure  8.  Experimentally-obtained  mean  plate  surface  temperature  shown  as  a  function  of 
time 


0i - , - 1 - 1 - 1 - i - 

0  10  20  30  40  50  60 

Time  (minutes) 


- 50-1 

- 50-2 

- 65-1 

- 65-2 

75-1 

75-2 


Figure  9.  Experimentally-obtained  maximum  plate  surface  temperature  shown  as  a 
function  of  time 
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In  an  effort  to  determine  the  influence  of  thermal  boundary  conditions  on  the  temperature 
increase  of  the  mock  energetic  plate  samples  when  subjected  to  harmonic  excitation,  a  thermally- 
insulated  box  was  constructed  to  roughly  simulate  an  insulated  environment.  The  structural 
integrity  of  the  box  was  maintained  by  acrylic  panes  with  a  3.8  cm  clearance  relative  to  the  top 
surface  of  the  plate,  and  additional  thermal  insulation  was  provided  by  1.9  cm  thick  foam 
insulation,  which  was  attached  on  all  sides.  The  box  attached  to  the  outside  of  the  clamping 
fixture  of  the  plate  sample  and  moved  with  the  shaker  head  during  the  course  of  experimentation. 

Due  to  the  nontransparent  nature  of  the  insulation  box,  the  temperature  profile  as  a  function  of 
time  could  not  be  obtained  with  the  infrared  thermal  camera  as  before.  As  a  result,  the  thermal 
camera  was  used  to  collect  a  snapshot  of  data  at  the  beginning  and  end  of  the  60  min  time 
interval  and  the  mean  and  maximum  temperature  difference  between  these  two  instances  was 
calculated.  A  similar  calculation  was  performed  on  the  ambient  temperature  profiles  presented 
above  to  yield  the  values  detailed  in  Table  3. 

Table  3.  A  comparison  of  the  experimentally-obtained  plate  surface  temperature  increases 
under  insulated  and  ambient  thermal  boundary  conditions  in  response  to  a  2  g 
excitation  near  the  first  resonant  frequency  for  all  of  the  plates 


Solids  Loading 
Sample 

Mean  Temperature  Increase 

Maximum  Temperature  Increase 

Insulated 

Ambient 

Insulated 

Ambient 

50%  -  1 

1.8  °C 

0.6  °C 

2.5  °C 

1.4  °C 

50%  -  2 

1.6  °C 

O 

oo 

o 

n 

2.1  °C 

1.7  °C 

65%  -  1 

1.2  °C 

0.3  °C 

1.5  °C 

1.0  °c 

65%  -  2 

1.0  °C 

0.1  °C 

1.6  °C 

1.0  °c 

75%  -  1 

2.1  °C 

0.3  °C 

2.5  °C 

1.0  °c 

75%  -  2 

1.9  °C 

0.2  °C 

2.3  °C 

0.7  °C 

As  evident  from  the  values  provided  above,  the  insulated  boundary  condition  does  result  in  a 
greater  temperature  increase  in  all  of  the  cases  as  expected.  However,  the  observed  trends  appear 
to  be  inconsistent  with  simple  bulk-scale  heat  transfer  models.  Current  efforts  are  focused  on 
determining  whether  this  is  due  to  experimental  error  (i.e.,  measuring  small  temperature  changes 
via  an  infrared  camera),  or  directly  attributable  to  the  particulate  composite  nature  of  the 
material. 

In  order  to  determine  an  analytical  (albeit  conservative,  due  to  the  nature  of  the  material)  upper 
bound  for  the  expected  temperature  increase  under  ambient  and  insulated  conditions,  a  heat 
generation  simulation  was  conducted  in  COMSOL  using  techniques  detailed  in  prior  work  [4]. 
The  mean  and  maximum  temperature  increases  for  the  simulated  plate  are  shown  in  Figures  10 
and  11  for  ambient  and  insulated  boundary  conditions.  Simulations  were  performed  on  both  50% 
and  75%  solids  loading  plate  samples  for  the  sake  of  comparison.  It  should  be  noted  that  though 
the  temperature  increases  reported  here  are  comparatively  small  (<10  °C),  they  can  be 
exacerbated  by  local  stress  concentrations  (a  focus  of  Task  Order  0002).  In  addition,  recent 
work  has  demonstrated  that  this  low-order  heating  (attributable  to  viscoelastic  loss  mechanisms) 
may  be  a  strong  indicator  of  a  sample’s  probability  of  achieving  thermal  run-away  and 
deflagration  [27] 
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Figure  10.  Simulated  mean  plate  surface  temperature  shown  as  a  function  of  time 


Figure  11.  Simulated  maximum  plate  surface  temperature  shown  as  a  function  of  time 


Experiments  identical  to  the  ambient  thermal  boundary  condition  tests  described  above  were 
performed  on  plates  of  varying  solids  loading  and  additive  content  created  as  part  of  the  PBXN- 
109  sample  set.  Note  that  the  data  presented  below  does  not  constitute  the  entire  data  set  as  there 
are  an  additional  six  plates  to  be  fabricated  and  tested  in  coming  months. 

In  order  to  minimize  the  effect  of  aging,  all  of  the  following  plate  samples  were  tested  at  an  age 
of  43  days.  The  vibrometer-recorded  HI  frequency  response  estimators  for  a  representative  plate 
in  ambient  conditions  in  response  to  three  levels  of  excitation  are  presented  in  Figures  12  and  13 
for  the  85%  solids  loading  plates  with  0%  and  30%  additive  content,  respectively. 
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Figure  12.  Experimental  HI  mechanical  frequency  response  estimator  obtained  for  an  85% 
solids  loading  -  0%  additive  content  plate  at  three  levels  of  excitation.  The  blue, 
green,  and  red  curves  depict  responses  at  1, 1.86,  and  2.44  g  RMS,  respectively. 
Solid  lines  represent  data  from  the  geometric  center,  and  dashed  lines  represent 
data  from  the  offset  point. 
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Figure  13.  Experimental  HI  mechanical  frequency  response  estimator  obtained  for  an  85% 
solids  loading  -  30%  additive  content  plate  at  three  levels  of  excitation.  The 
blue,  green,  and  red  curves  depict  responses  at  1, 1.86,  and  2.44  g  RMS, 
respectively.  Solid  lines  represent  data  from  the  geometric  center,  and  dashed 
lines  represent  data  from  the  offset  point. 


As  before,  the  corresponding  thermal  response  for  each  of  the  plate  samples  was  recorded  using 
the  FLIR  infrared  camera.  This  experiment  was  performed  over  a  60  min  time  interval  under  a 
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2  g  sinusoidal  excitation  near  the  first  resonant  frequency  of  each  respective  plate.  The  plate 
surface  temperatures  are  plotted  in  Figures  14  and  15.  The  legend  indicates  the  solids  loading 
content,  followed  by  the  additive  content  and  the  sample  number.  The  colored  envelope 
surrounding  each  curve  indicates  one  standard  deviation  for  each  trial.  Of  particular  note  here  is 
the  non-trivial  increase  in  the  observed  temperature  rise  in  the  presence  of  the  metallic  additive 
(aluminum). 


Figure  14.  Experimentally-obtained  mean  plate  surface  temperature  shown  as  a  function  of 
time 
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Figure  15.  Experimentally-obtained  maximum  plate  surface  temperature  shown  as  a 
function  of  time 
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4.2  Dissipative  Modeling  and  Material  Property  Parameter  Estimation 

As  alluded  to  above,  a  summary  of  the  random  base  excitations  testing  conducted  on  the  pure 
binder  and  high  solids  (85%  solids  where  100%  of  solids)  loading  sample  are  given  in  Figure  A1 
in  the  Appendix. 


A  sample  result  set  from  the  testing  on  a  92-days-old,  76.2  mm  (3  in)  tall/diameter  sample  of  the 
pure-binder  PBXN-109  is  shown  in  Figure  16.  In  this  test,  the  base  excitation  was  2  g  and  the 
sweep  was  from  8  to  68  Hz  in  60  s.  This  sample,  which  is  much  younger  than  the  previously- 
tested  pure-HTPB  samples  (only  92  days  after  fabrication)  has  stronger  nonlinear  behavior  than 
was  observed  with  the  older,  pure-HTPB  samples  which  were  more  than  180  days  old  when 
tested.  A  spectrogram  of  the  acceleration  response  of  the  mass  is  shown  in  Figure  16(a).  There 
are  clearly  harmonics  present  in  the  response;  the  harmonic  at  twice  the  excitation  frequency  is 
prominent.  A  harmonic  extraction  process  was  developed  to  examine  the  envelope  of  the 
response  at  the  driving  frequency,  and  the  results  of  applying  this  process  are  shown  in  Figure 
16(b).  Again,  some  evidence  of  quadratic  nonlinearity  is  seen;  there  is  an  increase  in  response 
from  the  general  trend  lines  at  half  the  resonance  frequency.  This  sample  is  much  more  lightly 
damped  and  softer  than  the  older  pure-HTPB  samples.  Estimates  of  the  linear  damping 
coefficient  and  Young’s  modulus  were  an  order  of  magnitude  lower  than  those  for  previously- 
tested,  older  pure-HTPB  samples. 
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Figure  16.  (a)  A  spectrogram  of  the  response  of  the  mass-material  system  undergoing  swept 
sine  excitation,  (b)  The  response  at  the  excitation  frequency.  The  sample  is  a 
76.2  mm  (3  in)  sample  of  HTPB  with  0%  solids  loading  excited  at  2  g. 
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As  an  example,  the  results  of  fitting  models  to  data  from  3  swept  sine  tests  on  a  76.2  mm  (3  in) 
cylinder  of  HTPB  50%  are  shown  in  Figure  17.  In  these  tests,  the  sample  weight  was  0.4  kg  and 
the  mass  loading  was  0.98  kg.  The  models  contain  either  a  linear  viscous  or  nonlinear  hysteretic 
damping  term,  an  A1'1  order  nonlinear  stiffness  term,  and  a  hereditary  viscoelastic  term,  which  is 
the  convolution  of  a  forcing  function  and  a  sum  of  exponentials  kernel.  Various  types  of 
viscoelastic  models  are  considered,  including  those  with  different  kernel  orders  and  different 
forcing  functions.  The  forcing  is  either:  (a)  a  function  of  the  relative  velocity  of  the  mass 
(directly  related  to  strain  rate);  (b)  a  function  of  the  nonlinear  stiffness  term  (directly  related  to  a 
nonlinear  function  of  strain);  or  (c)  a  function  of  the  hysteretic  damping  term.  The  model 
structures  that  were  used  to  produce  the  results  in  Figure  17  are  given  in  Table  4.  The  data  used 
in  the  estimation  of  the  model  parameters  were  from  three  swept  sine  tests  with  excitation  levels 
at  5  g,  7.5  g  and  10  g. 

Also  shown  in  Figure  17  are  the  envelopes  of  the  acceleration  of  the  mass  relative  to  the  base 
when  the  mass-material  system  is  undergoing  a  swept  sine  excitation  at  5  g.  Figure  17(b)  is  the 
part  of  Figure  17(a)  near  the  primary  resonance  of  the  system  (230  Hz  -  270  Hz).  The 
experimental  relative  acceleration  signal  envelope  (blue)  and  the  predicted  envelopes  from  the 
estimated  models:  Models  1  (green),  4  (red),  6  (black)  and  7  (purple)  are  shown.  All  of  the 
models  with  the  higher-order  viscoelastic  kernels  perform  better  than  Model  1  close  to 
resonance,  but  Model  1  performs  better  well  below  resonance.  In  this  region  (160-180  Hz) 
Models  4  and  6  over-predict  and  Model  7  under-predicts  the  relative  response  levels.  Above 
resonance,  Model  6  under-predicts  the  experimental  response  envelope,  whereas  the  responses  of 
Models  4  and  7  are  close  to  the  experimental  response  envelope. 


Table  4.  Model  numbers  and  corresponding  model  structure 


Model 

Number 

Stiffness 

Polynomial  Order 

(AO 

Damping 

iji 

l  Term 

a 

l 

Viscoelastic 
Kernel  Order 

(M) 

Viscoelastic 

Forcing 

1 

5 

Linear  X  =  0 

1 

Type  (a) 

4 

5 

Linear  X  =  0 

20 

Type  (a) 

6 

5 

Nonlinear  X  =  1.8 

20 

Type  (c) 

7 

5 

Nonlinear  X  =  1.8 

20 

Type  (a) 
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Figure  17.  (a)  The  experimentally-obtained  relative  acceleration  signal  envelope  (blue)  and 
the  predicted  envelopes  from  the  estimated  models:  Model  1  (green),  4  (red),  6 
(black)  and  7  (purple),  (b)  The  information  in  sub-figure  (a)  expanded  close  to 
the  resonance  of  the  system.  The  sample  here  is  HTPB  with  50%  solids  loading 
excited  at  5  g;  estimation  data:  3  swept  sine  tests  at  5  g,  7.5  g  and  10  g  base 
excitation  levels. 


4.3  Endochronic  Plasticity  Model  Results 

The  3D  endochronic  numerical  implementation  (detailed  in  depth  in  the  Appendix)  was  used  to 
model  the  loading  behavior  of  OFHC  copper,  and  the  model  results  were  compared  with 
experimental  results  reported  by  Lamba  and  Sidebottom  [28,29].  The  kernel  and  hardening 
parameters  were  calibrated  using  the  uni-axial,  constant  strain  amplitude,  cyclic  loading 
experimental  data  reported  in  [28].  Figure  18(a)  shows  the  prediction  of  stress-strain  response  of 
OFHC  copper  under  a  uni-axial  constant  strain  amplitude  cyclic  loading,  while  Figure  18(b) 
shows  the  prediction  of  the  shear  stress-strain  response  of  OFHC  copper  subjected  to  a  biaxial 
loading,  both  compared  with  experimental  results.  The  prediction  of  uniaxial  cyclic  loading  in 
Figure  18(a)  is  extremely  accurate,  thus  verifying  our  calibration  procedure  and  numerical 
implementation,  since  the  model  parameters  were  calibrated  directly  from  uniaxial  cyclic  loading 
data.  The  prediction  of  shear  stress-strain  response  in  Figure  18(b),  though  not  completely 
accurate,  gives  a  good  representation  of  the  experimental  measurements,  thus  validating  our 
calibration  procedure  and  numerical  implementation. 
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Strain  [%] 

(a) 


(b) 

Figure  18.  (a)  Uniaxial  constant  strain  amplitude  cyclic  hardening  response  of  OFHC 
copper  [endochronic  model  (solid  line)  vs  experimental  data  (symbols)],  (b) 
Prediction  of  the  shear  stress  versus  shear  strain  of  biaxial  non-proportional 
loading  of  the  OFHC  copper  [endochronic  model  (solid  line)  vs  experimental 
data  (dashed  line)]. 
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4.4  Crystal-Binder  Interface  Results 

The  elastic  constants  of  the  glass  bead  used  in  the  aforementioned  phase  field  simulations  are 
listed  in  Table  5.  To  choose  the  elastic  constants  for  Sylgard,  quasistatic  loading  simulations 
were  performed  and  compared  to  the  first  deformation  stages  of  the  experimental  force 
displacement  curve  shown  Figure  19.  The  best  fit  occurs  for  E  =  15  MPa. 

Table  5.  The  elastic  constants  for  Sylgard  and  glass  used  in  the  simulations 


■ 

Sylgard  184 

Glass  bead 

E 

75  MPa 

70  GPa  [30] 

V 

0.44  [31] 

0.3 

Figure  19.  Force  displacement  curves  used  for  the  calibration  of  the  elastic  constant  of 
Sylgard.  The  experiments  were  performed  by  the  Son  and  Chen  groups  at 
Purdue  University. 


Figure  20  shows  the  simulated  damage  field.  A  larger  amount  of  damage  is  observed  in  xz  plane 
due  to  the  asymmetry  of  the  dimensions  of  the  sample.  Figure  21  shows  the  delaminated  region 
in  the  experiments  and  in  the  simulations  for  an  applied  displacement  of  624  pm. 


Figure  20.  A  contour  plot  of  the  damage  field  for  an  applied  displacement  of  800  pm 
displacement  on  top  of  the  domain 
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Figure  21.  (a)  Experiment  and  (b)  simulation  of  debonding  for  a  624  pm  displacement 
applied  vertically 


5.  CONCLUSIONS 

The  research  described  herein  has  advanced  the  understanding  of  the  near-resonant  response  of 
energetic  materials  to  periodic  loading  through  a  joint  analytical,  numerical,  and  experimental 
investigation.  The  refined  modeling  tools  developed  with  this  improvement  in  basic 
understanding  are  currently  being  distilled  and  transferred  to  the  Air  Force  Research  Laboratory 
for  Department  of  Defense  use.  Despite  what  are  believed  to  be  positive  advancements, 
considerable  work  remains.  In  particular,  additional  research  focus  is  warranted  for  systems  with 
macroscopic  geometric  discontinuities/stress  concentrations  (e.g.,  cracks,  steps  in  geometry,  and 
intentional  holes),  alternate  sample  geometries,  and  less  traditional  binder  systems  (e.g.,  those 
with  a  high  degree  of  compliance  or  the  converse)  than  those  considered  here.  In  addition,  there 
is  a  need  for  continued  development  of  the  endochronic  material  model,  preliminary 
experimental  validation  of  this  model,  and  an  exploration  of  how  this  model  can  be 
mathematically  coupled  with  viscoelastic  models.  Finally,  efforts  should  be  made  to  distill 
lessons  learned  from  crystal-scale  computational  models  into  structural-scale  models  via 
"effective  material  properties"  or  an  alternate  multi-scale  approach. 
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A3.  Test  Schedule  and  Test  Reports 


The  information  in  this  appendix  relates  to  material  in  Sections  3.3  and  4.2  of  this  report.  It 
contains  the  base  excitation  material  testing  schedule  and  an  example  of  a  testing  summary 
report  for  a  random  base  excitation  test.  These  types  of  reports  are  automatically  generated  when 
a  test  is  completed  and  they  serve  as  a  record  of  the  test.  The  format  for  a  swept  sine  test 
summary  is  under  revision  because  in  recent  tests  responses  have  included  much  stronger 
harmonic  content  than  previously  seen  (Note:  we  are  now  testing  much  younger  samples  than 
previously  utilized,  as  well  as  a  slightly  different  composition)  and  we  would  like  to  capture  that 
information  in  the  test  report. 

The  geometry,  composition,  and  cure  date  of  the  samples  being  tested  along  with  the  testing 
progress  is  shown  in  Table  Al.  The  testing  matrix  is  divided  into  harmonic  and  random  base 
excitation  testing.  For  the  76.2  mm  (3  in.)  tall/diameter  high  solids  loading  sample  (85%  solids 
loading  where  70%  of  solids  is  sucrose  and  the  remaining  30%  is  the  aluminum  additive).  When 
a  series  of  10  scheduled  tests  is  terminated  due  to  the  fracturing  of  a  sample,  that  is  denoted  by  F. 
After  the  initial  10  scheduled  harmonic  and  random  tests  per  samples  are  completed,  we  will 
continue  tests  based  on  the  experimental  results  to  further  understand  the  dynamic  behavior  of 
the  samples. 

A  summary  report  of  a  random  base  excitation  test  is  shown  in  Figure  Al.  The  report  contains 
the  geometry  of  the  sample,  testing  conditions,  base  excitation,  and  response  of  the  signal,  as 
well  as  the  power  spectral  densities  and  the  magnitude  and  phase  of  the  estimated  frequency 
response  function  relating  the  base  acceleration  to  the  relative  acceleration  of  the  mass.  The 
estimated  coherence  is  also  shown.  The  frequency  response  estimates  for  the  random  excitation 
tests  are  estimated  by  using  the  HI  estimator  and  the  cross  and  power  spectral  density  estimates 
are  estimated  by  using  segment  averaging.  In  the  segment  averaging  a  Hann  window  was  used, 
segments  overlapped  by  50%  and  they  were  8192  points  long,  giving,  with  a  sample  rate  of  6000 
samples/s,  and  a  frequency  resolution  of  just  under  1  Hz.  Results  from  146  segments  are 
averaged  in  the  estimation. 
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Table  Al.  Summary  of  the  experimental  testing  status  of  the  PBXN-109  samples.  Here,  F 
denotes  Fracture,  and  Day  0  corresponds  to  the  sample  cure  date.  Each  sample 
has  its  own  number. 


Random  Base  Excitation  Harmonic  Base  Excitation 


Sample  ID 

0.015  0.025 

(g2/Hz)  (g2/Hz) 

(Day  #,  Day  #, 

0.050 

(g2/Hz) 

,...) 

Completed/ 

Total 

Planned 

F 

lg 

1 

2g  3g 

[Day  #,  Day  #,...) 

4g 

Completed/ 

Total 

Planned 

F 

P00_Q000_Q00_C0 1 
HTPB  100% 
Solids  0% 

Diam:  1"  Height:  1" 
Cure  Date:  2/29/2016 

77,86f 

2/10  F 

77 

77 

2/10  F 

P00_Q000_Q00_C02 
HTPB  100% 
Solids  0% 

Diam:  1"  Height:  1" 
Cure  Date:  2/29/2016 

86.102 

86,102 

4/10  F 

109 

109F 

2/10  F 

P00_Q000_Q00_C03 
HTPB  100% 
Solids  0% 

Diam:  3"  Height:  3" 
Cure  Date:  3/15/2016 

84,123 

84,123 

4/10 

85,92,99 

92 

85 

5/10 

P00_Q000_Q00_C04 
HTPB  100% 
Solids  0% 

Diam:  3"  Height:  3" 
Cure  Date:  3/28/2016 

78 

64,78 

78 

4/10 

79,86 

79,86 

79 

5/10 

P00_Q000_Q00_C05 
HTPB  100% 
Solids  0% 

Diam:  1"  Height:  1" 
Cure  Date:  6/20/2016 

4,25 

4,25 

4/10 

0/10 

P85_S070_A30_C03 
HTPB  15% 

Solids  85% 

Diam:  3"  Height:  3" 
Cure  Date:  6/13/2016 

15,29,36 

15,29,36 

6/10 

0/10 

P85_S  1 00_A00_C03 
HTPB  15% 

Solids  85% 

Diam:  3"  Height:  3" 
Cure  Date:  7/5/2016 

15 

15 

2/10 

0/10 

Note:  Sample  IDs  are  defined  by:  A##_B###_C##_D##.  A##  specifies  the  binder  formulation 
and  overall  solids  loading  (%  of  solid  in  overall  formulation).  For  example,  P  designates  that  the 
binder  is  HTPB  prepared  as  in  a  PBXN-109  propellant.  B###  specifies  the  type  of  mock 
energetic  crystal  and  percentage  of  solids  that  are  mock  energetic  (i.e.,  S100  means  all  solids  are 
sucrose).  C##  gives  the  type  of  additive  material  and  percentage  of  solids  that  are  additive. 
Here,  A  designates  Aluminum.  D##  gives  the  geometry  and  sample  identification  number. 
Here,  C  designates  a  cylinder.  The  Q  designation  followed  by  zeroes  in  any  place  essentially 
indicates  that  the  metric  is  not  applicable. 
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July  20, 2016  HTPB  85%  PBXN-109  C03 
Date  Curedjuly  5,  2016  Age:  15  days 
Sample  ID:  P85  S100  Q00  C03 


Diameter:  3  in 


Height:  3  in 
MaSSSample'  °  °  k'» 

MaSSMidplatc+‘  2  645 
Wavetek:  1500Hz 
Fs:  6000  Hz  T:  100  s  NFFT:  8192 
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Figure  Al.  A  summary  of  a  random  base  excitation  test  lasting  for  100  s  from  10  Hz  to 
1000  Hz  at  0.050  g2/Hz  for  sample  P85  S 100  Q00  C03 
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Coherence 


A4.  Modeling  and  System  Identification 

An  overview  of  the  mass-material  system  model  structure  and  the  techniques  used  to  estimate  the 
model  parameters  are  given  in  this  appendix.  The  mass-material  system  that  is  being  used  in  the 
experimental  tests  is  modeled  as  a  mass-spring-damper  with  a  viscoelastic  element.  Motion  is 
constrained  to  be  in  the  vertical  direction,  and  the  base  excitation  to  the  system  is  provided  by  a 
shaker.  The  experimental  apparatus  and  a  graphical  representation  of  the  system  model  are 
shown  in  Figures  A2. 


Figure  A2.  Experimental  test  setup  with  surrogate  sample  and  a  graphical  representation 
of  the  mass-material  system 


A4.1.  Model  Structure 

With  reference  to  Figure  A2,  the  equation  of  motion  of  the  mass  relative  to  the  base  is: 

mzj  +  c'^z^z-l)  +  k'iz-^  +  u1(/(z1,  Zi))  =  —my  —  mg  (Al) 

where  m  is  the  mass,  which  is  one-third  of  the  mass  of  the  material  plus  the  mass  of  the  center 
plate  and  any  masses  attached  to  it .  zx  —  x  —  y  is  the  relative  displacement  of  the  attached  mass, 
k' (zx)  is  the  stiffness  term,  c\z1,z1 )  is  the  damping  term,  v(f(z1,z1)')  is  the  viscoelastic  term, 
and  y  is  the  base  acceleration.  This  equation  can  be  expanded  about  the  static  settling  point  z0 
and  expressed  in  terms  of  the  motion  about  this  settling  point:  zx  =  z0  +  z.  Noting  also  that 
c'(0, z0)  +  /c'(zo)  +  t?(/(z 0,0))  =  —  mg  results  in  an  equation  about  the  settling  point  of  the 
form: 


mz  +  c(z,z)  +  k{z )  +  ^(/(z,  z))  =  —my  (A2) 

The  parameters  in  this  equation  are  functions  of  the  static  settling  point.  Thus,  to  extract  material 
parameters  from  estimates  of  the  parameters  of  this  equation  requires  knowledge  of:  (a)  the  static 
settling  point,  (b)  the  relationships  between  the  terms  in  Equation  Al  and  the  terms  in  Equation 
A2,  and  (c)  the  relationship  between  the  parameters  in  Equation  Al  and  the  parameters  of  the 
corresponding  stress-strain  form  of  that  equation.  In  this  appendix,  the  focus  is  on  describing 
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techniques  to  estimate  the  parameters  of  Equation  A2,  but  it  is  noted  that  several  additional  steps 
are  required  to  translate  those  estimates  into  material  parameters. 

Many  different  forms  of  the  damping,  stiffness,  and  viscoelastic  terms  are  considered  in  the 
modeling  process.  The  terms  may  depend  on  the  relative  velocity  (z),  the  relative  displacement 
(z),  or  a  combination  of  both.  From  linear  viscous  damping  (cz)  and  linear  stiffness  ( kz ),  to 
more  complex  expressions,  such  as  so-called  hysteretic  damping  terms  of  the  form  (cz\z\A). 
When  considering  nonlinear  stiffness  models,  an  Nth  order  polynomial  expression  in  z 
(Zi=i  k-iZ1)  is  used  because  it  is  linear  in  parameters.  Use  of  Ogden's  hyperelastic  model  has  been 
explored  previously  [Al,  A2]  and  Paripovic  and  Widdle  both  showed  that  an  Ogden  model  can 
be  fitted  to  the  estimated  polynomial  stiffness  term  to  extract  hyperelastic  material  parameters. 
For  the  viscoelastic  term,  a  hereditary  model  was  chosen.  The  general  form  of  this  model  is  a 
convolution  of  a  hereditary  (sum  of  exponentials)  kernel  and  a  forcing  function: 

v(t)  =  f(z,z )  *  g(t)  =  f  f(z(z),z(z))g(t  -  z)dz  (A3) 

J  —  CO 

where  g{t)  is  the  relaxation  kernel,  /(z,  z)  is  an  arbitrary  driving  function  that  may  depend  on 
the  relative  displacement,  relative  velocity  or  both.  The  relaxation  kernel,  g{t)  is  of  the  form: 

M 

g(t)=YjPi^““  (A4) 

i= 1 

where  M  is  the  order  of  the  viscoelastic  model,  and  cq  and  fa  are  the  viscoelastic  model 
parameters.  Combining  Equations  A3  and  A4  gives: 

t  M 

v(t)  =  f(z,z)  *  g(t)  =  f  f(z(z),z(z))\'fae~ai(t~^dz  (A5) 

J~°°  i=i 

When  modeling  the  dynamic  behavior  of  the  HTPB  samples,  first,  the  linear  case  was  considered 
with  no  viscoelastic  term  and  then  terms  were  systematically  added  and/or  adapted  to  improve 
the  accuracy  of  the  model  predictions.  The  various  model  structures  that  have  been  examined  so 
far  are  shown  in  Table  A2. 
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Table  A2.  All  of  the  stiffness,  damping  and  viscoelastic  term  combinations  considered  when 
fitting  the  mass  material  model  to  the  simulated  responses  and  experimental 
data  from  the  sweep  and  random  base  excitation  testing  on  HTPB  surrogate 
samples 


Model 


Stiffness 

Damping 

Viscoelastic  Term 

Kernel  Order 

k(z ) 

c(z) 

f(z,  z) 

M 

/qz 

cz 

- 

- 

V  ktzl 

cz 

- 

- 

■*— <i=l 

v5 

cz|z|A 

/  kiZl 

y  kizi 

cz 

z(t) 

1 

Z—>i= l 

y  kizi 

cz 

z(t) 

various 

Z—ti= 1 

y  ^z1 

cz 

y  ^z1 

1 

Z—>i= l 

Z—>i= l 

y  ^z1 

cz|z|A 

z|z|A 

various 

Z—>i= l 

y  ^z1 

cz|z|A 

z(t) 

various 

Z—ti.  =  1 

1 

2 

3 

4(a) 

4(b) 

5 

6 
7 


A4.2.  Model  Parameter  Estimation  for  Low  Order  Viscoelastic  Models 

In  [3]  Paripovic  focused  on  models  with  no  viscoelastic  term  (v  —  0)  or  with  a  viscoelastic 
kernel  of  order  1  (M  —  1).  The  estimation  technique  used  for  these  cases  was  continuous-time 
system  identification.  When  viscoelastic  kernels  are  higher  order,  the  method  described  in  this 
section  becomes  difficult  to  use,  and  so  an  iterative  parameter  estimation  technique  was 
developed  for  when  M  >  2;  that  technique  is  described  in  a  later  section.  The  iterative  method 
also  utilizes  continuous-time  estimation  to  estimate  the  stiffness  and  damping  model  parameters. 

A4.2.1.  Continuous-Time  System  Identification  Approach 

The  continuous-time  system  identification  methodology  is  introduced  with  its  application  to 
Model  2,  where  a  linear  viscous  damping  term,  a  5th  order  polynomial  stiffness,  and  v(t)  —  0  is 
considered.  The  equation  for  Model  2  is: 

mz  +  cz  +  kxzx  +  k2z2  +  k3z3  +  k4z4  +  k5z5  —  —my  (A6) 

Note  that  Equation  A6,  while  nonlinear  in  z,  is  linear  in  parameters  (/c,  and  c).  The  mass  m  is 
known,  it  is  a  sum  of  the  mass  of  the  midplate  (0.98  kg),  the  thin  aluminum  plate  that  is  glued  to 
the  sample,  one  third  of  the  test  specimens  mass,  and  any  additional  mass  fixed  to  the  midplate. 
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Moving  all  of  the  known  quantities  to  the  right  hand  side  and  using  the  time  histories  samples, 
Equation  A6  may  be  written  as  a  system  of  equations  of  the  form: 


— m(z(A)  -1-  y(A)) 

"  z(A)  z(A)  z2(A) 

.  z5(A)  ■ 

r  c 

ki 

— m(z(2A)  +  y(2A)) 

= 

z(2A)  z(2A)  z2(2A)  . 

.  z5(2A) 

k-2 

/e  3 

m(z(NA)  +y(NA))_ 

_z(AA)  z(AA)  z2(AA) 

••  z5(1VA)_ 

k4 

-k-5- 

This  can  be  expressed  as: 

q=[A]p  (A7) 

where  A  is  the  time  step  in  seconds.  Note  that  each  row  in  A  and  q  is  time  aligned  [A4].  Here,  A 
and  q  are  known  since  the  experimental  signals  are  measured  and  can  be  manipulated  and 
conditioned  appropriately  to  yield  z2,  z3,  etc.  p  is  the  vector  of  the  unknown  terms.  One 
advantage  of  this  method  is  that  the  time  domain  signals  may  be  equally  cropped  to  focus  on  a 
specific  region  of  the  response  and  avoid  unwanted  transient  or  noise  corrupted  portions  of  the 
signal.  Solving  the  normal  equations: 

{Arq}  =  [ATA]p  (A8) 

where  T  denotes  transpose,  yields  an  estimate  of  the  parameter  vector  p.  This  can  be  solved,  for 
example,  by  using  the  regress  command  in  MATLAB,  which  also  provides  confidence  intervals 
and  statistics  related  to  the  fit. 

A4.2.2.  Signal  Processing 

Before  applying  this  approach  to  the  measured  experimental  signals  the  signals  must  be 
appropriately  filtered  and  conditioned  to  reduce  ill-conditioning  of  the  matrices  which  can 
impact  the  accuracy  of  the  parameter  estimates. 

One  challenge  of  applying  this  approach  to  experimental  measurements  is  deriving  the  necessary 
velocity,  displacement,  and  other  higher  order  state  space  signals  from  the  measured  relative 
acceleration  (z).  One  must  be  careful  to  condition  the  signal  appropriately  to  avoid  high  levels  of 
noise  corruption  on  the  signals.  Also,  care  must  be  taken  to  make  sure  that  all  of  the  columns  of 
the  estimation  matrix  contain  time  histories  that  are  equally  filtered.  For  example,  when  deriving 
estimates  of  the  velocity  and  displacement  signals  from  the  measured  acceleration,  the  following 
process  is  used.  First,  the  measured  relative  acceleration  signal  is  passed  through  a  Butterworth 
high  pass  filter  with  a  cut-off  frequency  and  roll-off  rate  chosen  by  considering  the  excitation 
range  and  the  response  characteristics  of  the  system.  The  signal  is  then  integrated  by  using  a 
digital  filter  integrator  to  yield  an  estimate  of  the  relative  velocity.  This  two-step  process  is 
repeated  once  more  to  yield  an  estimate  of  the  relative  displacement  measurement.  To  ensure 
equal  signal  conditioning  of  the  time  histories  used  in  the  estimation  matrix,  the  relative 
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acceleration  and  base  acceleration  are  also  passed  through  the  same  high  pass  filter  twice  and  the 
estimated  relative  velocity  once  more. 

The  high  pass  filters  and  integration  introduce  unwanted  transients  at  the  start  and  ends  of  the 
signal.  These  parts  of  the  signal  are  removed.  Typically,  the  first  and  last  10%  of  the  signals  are 
cropped.  Further  cropping  of  the  signal  may  be  used  to  use  only  specific  regions  of  the  time 
history  that  may  be  of  interest  (e.g.,  near  resonance)  in  the  estimation.  The  newly  conditioned 
relative  acceleration,  velocity,  displacement,  forcing,  and  other  terms  are  now  ready  to  be  used  in 
the  continuous-time  system  identification  approach.  When  higher  order  models  are  being 
considered,  derivatives  of  the  acceleration  signals  are  needed.  In  this  case,  a  combination  of  low- 
pass  filters  to  reduce  high  frequency  noise  and  digital  filter  differentiators  are  used,  and  all  of  the 
signals  receive  the  same  amount  of  low-pass  filtering,  as  with  the  high-pass  filtering  described 
above. 

A4.2.3.  Summary  of  the  Continuous-Time  Process 

The  general  approach  of  applying  the  continuous-time  system  identification  to  estimate  the 
system  parameters  is  as  follows: 

1.  Based  on  the  physical  system,  derive  an  equation  of  motion  with  the  appropriate 
stiffness,  damping,  and  viscoelastic  terms  to  model  the  desired  dynamic  behavior. 
The  equation  should  be  a  function  of  measurable  or  derivable  states  (velocity, 
displacement,  acceleration,  or  higher  derivatives)  and  the  model  parameters.  Ideally, 
it  will  be  linear  in  parameters. 

2.  Derive  the  needed  state  space  signals  defined  in  the  equation  of  motion.  Be  sure  to 
condition  the  signals  appropriately  and  filter  all  signals  by  an  equal  amount.  Filters 
should  be  designed  to  reduce  noise  corruption  and  signals  should  be  truncated  to 
remove  transient  filter  responses. 

3.  For  the  linear  in  parameters  case,  re-arrange  the  equation  of  motion  so  that  all  terms 
that  are  known  are  on  the  left  hand  side  and  all  terms  that  are  functions  of  the 
parameters  are  on  the  right  hand  side.  Construct  the  vector  q  and  the  matrix  A  using 
the  conditioned  signals. 

4.  Depending  on  the  structure  of  A,  the  matrix  may  become  ill-conditioned.  This  results 
in  MATLAB  using  a  pseudo-inverse  method  to  solve  the  system  of  equations  and 
setting  small  singular  values  to  zero,  but  these  small  terms  may  have  important 
physical  significance.  To  try  to  overcome  this,  it  is  recommended  that  the  individual 
columns  of  A  are  divided  by  their  respective  standard  deviation.  The  division  will 
normalize  the  matrix  such  that  all  of  the  columns  are  of  similar  order  and  usually  the 
estimation  accuracy  improves.  Note  that  if  the  matrix  A  is  normalized,  the  calculated 
vector  p  needs  to  be  re-normalized  to  yield  the  correct  stiffness  and  damping 
estimates. 

5.  By  using  a  linear  regression  program  (for  example,  regress  in  MATLAB)  or 
something  similar,  solve  the  normal  Equations  A7  to  yield  estimates  of  p.  When  the 
parameters  in  p  are  functions  of  the  physical  parameters,  it  will  be  necessary  to  define 
a  process  to  take  the  elements  of  p  to  derive  estimates  of  the  system  parameters  of 
interest,  for  example,  (klt  k2,  —  ,c ). 
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6.  Use  the  parameter  estimates  to  calculate  the  system  response  by  using,  for  example,  a 
Runge  Kutta  algorithm  with  an  appropriate  step  size.  Compare  the  estimated  response 
with  the  original  response  as  an  indication  of  the  accuracy  of  the  model.  The  difficult 
part  is  then  to  determine  what  model  structure  adjustments  are  needed  to  improve  the 
predictions. 

A4.2.4.  Inclusion  of  the  Viscoelastic  Terms  in  Continuous-Time  Estimation 


One  approach  when  v(t)  is  not  zero  is  to  reformulate  the  integro-differential  equation  into  a 
system  of  coupled  ordinary  differential  equations.  The  continuous-time  system  identification 
approach  can  then  be  applied  to  estimate  the  system  parameters.  To  reformulate  the  viscoelastic 
term,  take  the  Laplace  transform  of  Equation  A5 

M 

v{t)=Yjfit)*Pie~ait  (A9) 

i= 1 


and  note  that  the  Laplace  transform  of  a  convolution  (in  time)  is  a  multiplication  in  the  s 
domain.  Thus,  the  Laplace  transform  of  Equation  A9  is  of  the  form: 


M 


L(s)=^/(s) 


Pi 


1=1 


s  +  CLi 


=  F(S ) 


A%-iSm  1  +  -  +  M0 
sM  +yM_1sM~1  +  -  +  Yo 


(A10) 


Multiplying  through  by  the  denominator  of  Equation  A 10,  rearranging  and  taking  the  inverse 
Laplace  transform  yields  the  ordinary  differential  equation: 


dMv  dM  1v  dM  1/(t) 

+  Ym- i  +  •"  +  YoV  =  Hm-i  dtM- 1  +  •"  +  Mo/(0 


(All) 


where  [it  and  yt  are  combinations  of  cqand  /^. 


Prior  to  being  able  to  apply  the  continuous-time  system  identification  method  to  estimate  kt,  c, 
lii  and  Yi  (and  then  cqand  /?j),  Equation  A2  is  rearranged  to  express  v(t)  as  a  function  of  terms 
that  are  themselves  a  function  of  z  and  its  derivatives.  It  is  then  differentiated  M  times  and 
substituted  back  to  build  an  (M  +  2)th  differential  equation.  Solving  Equation  A2  for  v(t),  and 
differentiating  yields: 


and 


v(t)  —  —my  —  mz  —  /c(z)  —  c(z), 
i>(t)  =  -my  -mz-jt  [/c(z)]  -  ^  [c(z)]. 


(A  12) 


(A13) 


etc. 
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Substituting  Equation  A12,  Equation  A13,  etc.  into  Equation  A4  results  in  the  (M  +  2)th 
differential  equation  in  z.  For  example,  for  Model  4a  (with  M  —  1  and  N  —  1),  the  following 
system  of  equations  may  be  constructed: 


‘  — m(z(A)  +  y(A))  ‘ 

"  z(A)  z(A)  z(A)  m(y(A)) 

'  c  +  TTLCCi ' 

— m(z(2A)  +  y  (2A)) 

z(2A)  z(2A)  z(2A)  m(y(2A)) 

k  +  coq 

— m(z(3A)  -I-  y  (3A)) 

— 

z(3A)  z(3A)  z(3A)  m(y(3A)) 

ka±  +  /?x 

m(z(NA)  +  y(NA))_ 

_z(AA)  z(AA)  z(NA)  m(y(NA))_ 

.  _ 

This  is  a  linear  system  of  equations,  similar  to  Equation  A7  but  with  the  added  complexity  of 
having  to  transform  the  estimated  parameters  to  estimates  of  c,  k,  cqand  /?,  .  Note  that,  in  general, 
up  to  (M  +  2)th  derivatives  of  the  state  space  signals  z  and  y  are  needed  to  construct  this  matrix, 
and  while  for  the  M  —  1  the  structure  remains  simple,  increasing  M  increases  the  number  of 
terms  significantly  and  there  are  also  signal  and  matrix  conditioning  problems  when  estimating 
high  order  derivatives.  To  avoid  these  problems,  an  iterative  estimation  approach  was  developed. 
This  is  described  in  Section  A4.3. 

A4.3.  Estimation  of  Higher  Order  Viscoelastic  Models 

The  previously  described  continuous-time  system  identification  approach  is  not  desirable  for 
higher  order  models  because  the  (M  +  2)th  derivatives  of  the  signals  are  needed  and  deriving 
these  from  measured  acceleration  signals  can  result  in  noise  corruption  of  the  signals  [A5,  A6]. 
To  improve  the  estimation  an  iterative  approach  was  developed.  The  parameter  estimation  is  split 
into  two  parts.  In  one  part,  the  stiffness  and  damping  coefficients  are  estimated  by  using  the 
continuous-system  identification  approach  described  in  the  previous  section,  and  in  the  other  part 
the  impulse  and/or  frequency  response  of  the  viscoelastic  kernel  is  modeled  by  using  Prony 
Series  estimation  techniques.  The  estimated  viscoelastic  term  v(t )  forms  the  bridge  between  the 
two  parts  of  the  estimation.  A  visual  illustration  of  the  iterative  process  is  shown  in  Figure  A3. 
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Figure  A3.  The  steps  in  the  iterative  system  identification  approach  used  for  parameter 
estimation.  The  black  outlined  boxes  illustrate  the  initialization  of  the 
estimation,  and  the  blue  outlined  boxes  are  the  core  of  the  iterative  estimation 
scheme. 


There  are  10  steps  in  this  estimation  process,  which  are  given  below.  Note  that  it  is  assumed  that 
the  base  excitation  frequency  range  and  amplitude  are  sufficient  to  provide  enough  information 
to  determine  the  system  parameters.  Note  that  it  is  straightforward  to  incorporate  signals  from 
different  tests  into  the  estimation  process,  which  is  helpful  when  trying  to  develop  a  single  model 
of  the  system  that  predicts  responses  to  different  levels  and  types  of  excitation.  In  reality,  the 
desired  excitations  will  evolve  as  more  is  understood  about  the  system  behavior,  model 
deficiencies  are  identified,  and  the  proposed  improved  models  and  their  identification  are 
explored  through  simulation  and  analysis.  Here,  just  the  parameter  estimation  process  is 
described  for  a  given  nonlinear  and  viscoelastic  model  structure. 

Step  1:  Import  measured  signals  from  experiments:  relative  acceleration  z  —  x  —  y  and  base 
acceleration  y.  Make  sure  to  subtract  any  offset  to  ensure  that  the  accelerations  are  centered 
about  zero. 


Step  2:  Integrate  the  relative  acceleration  signals  to  calculate  the  states  of  the  model,  z,  y  -*  z,z 
also,  depending  on  model  structure,  calculate  the  necessary  polynomial  expressions  of  the  signals 
(z2,  ...,  etc.). 

Step  3:  Assign  an  initial  value  to  the  viscoelastic  term  v(t).  Currently  we  are  assuming  the 
initial  value  of  v(t )  is  zero. 
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Step  4:  By  using  the  continuous-time  system  identification  approach,  estimate  the  stiffness  and 
damping  terms  (c;  kf,  i  =  1,  2,  ...N).  Note  that  if  high  order  polynomial  expressions  of  the 
stiffness  terms  are  considered,  it  will  be  particularly  important  to  pay  attention  to  the 
conditioning  of  the  system  matrix  of  Equations  A7,  prior  to  solving  using  the  normal  equations. 
This  is  accomplished  by  dividing  each  column  by  its  respective  standard  deviation.  Note,  to 
compensate  for  this  normalization,  the  estimates  in  the  parameter  vector  p  will  need  to  be 
multiplied  by  the  corresponding  standard  deviations  to  yield  estimates  of  the  model  parameters. 

Step  5:  Update  the  viscoelastic  time  history  estimate  by  subtracting  the  damping  and  stiffness 
terms  from  the  inertial  terms:  v(t )  =  —my  —  mz  —  c(z,z)  —  k(z). 

Step  6:  Calculate  the  viscoelastic  kernel  frequency  response  function,  G(fk ).  If  the  signals  are 
deterministic  G(fk)  —  DFT[v(nA)]/DFT[f(nA)],  or  if  the  signals  are  random  G(fk)  —  SfV/Sff. 
Note  that  DFT[— ]  denotes  a  Discrete  Fourier  Transform  and  Sab  denotes  the  estimated  cross 
spectral  density  between  signal  a  and  signal  b.  Segment  averaging  is  used  with  Hann  windows  at 
a  50%  overlap  to  estimate  the  cross  and  power  spectral  densities.  An  option  in  Step  7  is  to 
directly  fit  to  the  impulse  response  of  the  viscoelastic  kernel.  An  estimate  of  the  impulse 
response  is  produced  by  applying  an  Inverse  Discrete  Fourier  Transform  IDFT[— ]  to  the 
estimated  frequency  response  function  and  rearranging  to  place  any  acausal  behavior  (an 
outcome  of  all  the  filtering  and  signal  conditioning)  at  the  start  of  the  impulse  response. 

Step  7:  Fit  an  Mx\h  order  Prony  Series  model  (M1  »  M )  either  to  the  frequency  response 
function  or  to  the  impulse  response  function.  Calculate  the  energy  contributions  of  each  term 
over  a  specified  time  period  tx  to  t2  of  the  impulse  response.  This  time  is  specified  to  avoid  any 
acausal  behavior  at  the  start  and  to  not  include  too  much  of  the  region  where  the  impulse 
response  has  decayed  into  the  noise  floor.  In  the  estimation  process  it  is  possible  for  the 
viscoelastic  terms  to  be  either  real  or  complex.  Because  the  signal  is  real,  complex  terms  appear 
in  conjugate  pairs.  When  terms  are  complex,  the  energy  contribution  of  the  complex  conjugate 
pair  is  calculated.  By  energy,  we  mean  the  sum  of  the  squared  signal  from  tx  to  t2. 

Step  8:  Reduce  the  model  order  to  M  by  sorting  the  viscoelastic  terms  by  their  energy 
contribution  and  selecting  the  terms  with  the  highest  individual  energy  contributions.  Additional 
discussion  on  model  reduction  techniques  and  the  challenges  presented  are  included  in 
subsection  A4.3.1.  [Note  that  there  are  some  problems  where  two  terms  (or  four  complex  terms) 
nearly  cancel,  thus  their  net  contribution  is  small,  but  the  individual  contributions  are  large.  Also, 
complex  terms  are  sometimes  associated  with  frequencies  well  outside  the  excitation  region. 
Refinement  of  model  term  selection  is  a  part  of  on-going  research.] 

Step  9:  Check  the  convergence  criteria.  Currently  this  is  specified  to  be  a  less  than  0.1%  change 
in  the  parameter  estimates  for  the  stiffness  and  damping  terms,  and  a  less  than  1%  change  in  the 
estimated  viscoelastic  parameters.  These  criteria  are  preset  in  the  program.  Increasing  the 
tolerances  decreases  the  computational  time.  Note  that  if  the  convergence  criteria  are  too  tight, 
the  model  will  never  converge  and  the  parameter  estimates  will  continue  to  oscillate  between  two 
values  until  maximum  iteration  number  has  been  met. 
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Step  10:  If  convergence  is  not  achieved,  the  viscoelastic  term  estimate  is  updated  and  the  process 
returns  to  Step  4  and  repeats. 

A4.3.1.  Prony  Analysis  and  Model  Order  Reduction  Discussion 

There  are  many  methods  described  in  literature  that  can  be  used  to  reduce  the  model  order  of  an 
estimated  model,  specifically  with  application  to  estimated  Prony  Series  models.  One  main 
advantage  of  using  Prony  analysis  is  that  it  is  a  parametric  method  that  fits  a  sum  of 
exponentials,  which  is  of  the  same  structure  as  the  viscoelastic  hereditary  kernel.  When  applied 
to  a  measurement,  the  Prony  Series  models  both  the  signal  structure  and  the  noise  structure,  as 
well  as  any  other  artifacts  introduced  during  the  signal  processing,  and,  in  this  case,  impulse 
response  estimation.  Model  reduction  is  desired  because  while  a  higher  order  model  may  result 
in  an  excellent  fit  to  the  measurements  or  estimated  signals,  many  terms  may  have  nothing  to  do 
with  the  behavior  of  the  material.  A  simpler  model  that  models  the  viscoelastic  behavior  of  the 
material,  as  opposed  to  other  artifacts,  is  desirable  [A7,  A8].  Extracting  only  these  physical  terms 
can  be  challenging. 

Currently,  we  employ  only  individual  term  energy  contributions  as  the  criterion  for  developing 
an  Mth  order  viscoelastic  kernel  from  a  higher  order  Mjth  order  model.  The  energy  contribution 
terms  are  calculated  for  each  team  by  squaring  the  signal  and  integrating  over  the  same  time 
interval  for  the  system.  When  Prony  Series  terms  appear  in  complex  conjugate  pairs  the  energy 
contribution  of  the  pair  is  determined  by: 

rti 

E.C.=  I  (Ae~at  cos(2nft  —  <p))2dt  (A15) 

'ti 

and  for  a  real  term  the  energy  contribution  is  of  the  form, 

^2 

(A  e~at')2dt  (A16) 

1 

When  simulating  the  viscoelastic  term,  it  is  straightforward  to  choose  the  time  interval  for  the 
energy  contribution  integral;  however,  when  the  viscoelastic  term  profile  is  unknown  it  is  not  as 
trivial.  Ideally,  the  interval  is  chosen  such  that  the  initial  transient  behavior  is  ignored  and  the 
terms  have  decayed  enough  over  the  region  of  integration  so  that  the  amplitudes  are  in  the  noise 
floor.  A  very  long  interval  is  not  wanted  because  eventually  the  integral  will  calculate  energy 
contributions  of  the  noise  corruption.  The  nontrivial  part  is  that  the  time  interval  must  remain 
equal  for  all  the  contributing  terms  and  be  a  good  region  for  all  terms. 

As  an  example,  the  Prony  Series  analysis  is  applied  to  the  experimentally  determined  viscoelastic 
impulse  response  resulting  from  testing  on  a  pure-HTPB  76.2  mm  (3  in)  material  sample  with  a 
0.98  kg  mass  loading  (corresponding  to  m  =  1.12  kg).  This  system  was  excited  at  5  g.  After 
calculating  the  energy  contribution  of  each  term  of  the  =  250  Prony  Series  estimation,  the  top 
ten  contributing  terms  were  selected  and  used  to  recreate  the  experimental  impulse  response. 
This  is  shown  in  Figure  A4.  Even  with  10  terms,  with  the  reduced  order  model  it  is  possible  to 
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predict  the  first  few  oscillations  of  the  impulse  response  before  it  decays  to  the  noise  floor  in  the 
predicted  impulse  response. 
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Figure  A4.  The  impulse  response  of  a  pure-HTPB  sample  (blue)  excited  at  5  g  along  with 

the  predicted  Prony  approximation  for  (left)  =  250  (red)  and  (right)  the  top 
energy  contributing  terms,  M  —  10  (green) 


Another  approach  is  to  fit  directly  to  the  estimated  frequency  response  function.  There  are  parts 
of  the  frequency  response  function  that  are  heavily  corrupted  by  noise  and  parts  where  the 
estimation  is  reasonably  good  (in  the  region  of  excitation).  By  doing  the  estimation  in  the 
frequency  domain,  we  can  select  regions  of  the  frequency  response  where  the  information  is  not 
strongly  corrupted  by  noise. 

The  model  of  the  sampled  system  transfer  function  that  relates  the  sampled  forcing  function  to 
the  sampled  viscoelastic  function  is  determined  by  sampling  the  viscoelastic  impulse  response 
g(t )  and  taking  the  z-transform  of  g(nA).  It  is  of  the  form: 


G(z)  = 


b0  +  bxz  1  +  -  +  bM_xz  (M  1} 
1  +  a1z_1  +  — I-  aMz~M 


(A17) 


Substituting 


j2nf , 

z  —  e  fs 

where  /  is  the  frequency  and  fs  is  the  sampling  frequency,  gives  the  frequency  response  of  the 
model: 
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G(z)  = 


—j2nf  /  -;'27t/(M-1), 

b0  +  bxe  fs  +  •••  +  bjyj^iS  fs 


1  +  aqe 


-j2nf,  —j2nfM  / 

fs  -|-  •••  +  CLMe  fs 


(A18) 


Equation  A18  can  be  rearranged  and  evaluated  at  frequencies  /  =  fi,f2,  —  ,fN  and  written  in  the 
form: 


G(A  )' 
G(/2) 

-£(/*)- 


-j2nf1  , 

-GtJJe  jf*  . 

■■  1 

e  7 fs 

-jlnh/ 

-G(f2)e  't,  . 

■■  1 

~j2nf2  , 
e  7 fs 

-]2nfN/ 

-G(fN)e  'fs  ■ 

-j2nfN  , 

e  7 fs 

..  1 

e 

e 

e 


-j2n{M-l)f1/ 

'fs 

—j2n(M—l)f2  I 
'fs 


-j2n(M-l)fN  j  J 


aM 

b0 

bi 


lb 


M-  1J 


where  M  is  the  model  order,  fs  is  the  sampling  frequency  and  G(/i)  is  the  viscoelastic  frequency 
response  evaluated  at  fo  Hz.  This  can  be  written  more  simply  as: 


{G}  =  l X]{P } 


(A19) 


Note  that  estimation  errors  due  to  noise  and  un-modelled  dynamics  when  solving  Equation  A19 
can  sometimes  yield  complex  imaginary  coefficients.  To  prevent  this,  the  matrix  equation  is 
reformulated  to  force  the  imaginary  part  of  the  coefficients  equal  to  zero.  The  new  matrix 
equation  is  of  the  form: 


Re{G  (A)} 

\Re{X11] 

Re{X  i2}  . 

r  ai  1 

7m{G(/i)} 

lfn{X  n} 

MA12}  .  - 

a2 

Re{G(fN)} 

Im{G(fN)y 

Re{XN  A 
Im{XN1 } 

Re{XN2]  . 

Im{XN2 }  . 

2 

1- 

The  model  order  reduction  technique  applied  in  the  frequency  domain  is  less  sensitive  to  noise 
corruption  than  the  Prony  Series  applied  in  the  time  domain.  One  explanation  is  that  the  selected 
frequency  region  is  in  the  region  of  excitation  and  in  this  region  the  estimation  is  less  corrupted 
by  noise.  When  applying  Prony  Series  analysis  to  experimental  data  in  the  frequency  domain  it 
was  found  that  lower  order  viscoelastic  models  were  sufficient  to  model  the  system  behavior. 
The  results  of  fitting  a  Prony  Series  model  to  the  estimated  viscoelastic  frequency  response 
function  ( vFRF )  derived  from  testing  on  a  pure-HTPB  sample  excited  at  10  g  with  order  = 
40,  and  selecting  only  the  top  10  contributing  terms  for  varying  regions  of  estimation  is  shown  in 
Figure  A5. 
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Figure  A5.  The  viscoelastic  kernel  frequency  response  magnitudes  (\vFRF\)  derived  from 
experiments  on  pure-HTPB  excited  at  10  g.  The  gray  region  denotes  the  region 
of  strong  excitation  in  the  experiment,  red  is  \vFRF\  derived  as  part  of  the 
iterative  whole  model  estimation  process,  blue  is  the  region  used  in  the  vFRF 
estimation  process,  and  green  is  the  \vFRF\  predicted  using  the  estimated 
viscoelastic  kernel  model  when  fitting  M1  —  40  and  selecting  only  the  top  10 
terms.  The  regions  used  in  estimation:  (a)  50-110  Hz,  and  (b)  30-130  Hz. 


This  viscoelastic  kernel  estimation  and  model  order  reduction  has  been  incorporated  into  the 
model  estimation  process,  but  currently  the  user  still  has  to  specify  a  desired  model  order.  Future 
work  will  include  investigation  of  automatic  selection  of  model  order. 
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A5.  3D  Numerical  Implementation  of  Endochronic  Constitutive  Equations 

A5.1.  Endochronic  Constitutive  Equations  and  their  Numerical  Integration 

The  numerical  integration  of  the  endochronic  constitutive  equations  follows  the  method  reported 
by  Hsu  et  al.  [A9],  modified  to  include  rate-dependent  (viscoplastic)  effects. 

In  order  to  proceed  with  the  numerical  procedure,  it  is  imperative  to  assume  a  functional  form  of 
the  kernel  function  p(z).  According  to  Hsu  et  al.  [A9],  for  numerical  applications  the  kernel  can 
be  assumed  to  be  expressed  as  the  sum  of  a  group  of  n  exponentially  decaying  functions  of  the 
form: 


ii 

pOO  =  £ 


Cr  e  UrZ 


(A21) 


r= 1 


Substituting  the  above  approximation  for  p(z)  in  Equation  1,  we  get 


Zrz  ,  dev 

Y  J  Cre-ar(z-z  )-r-dz' 


r= 1 


dz' 


(A22) 


Let  the  loading  be  divided  into  m  steps,  and  the  following  notations  be  used: 

fe(  ):  value  of  a  variable  at  the  end  of  kth  loading  step  (k  =  0,1,2  ...  m). 
k(A):  increment  in  the  value  of  variable  caused  by  k,h  loading  step,  i.e.,  fe()  : 

Now,  the  integral  in  Equation  A22  can  be  expressed  approximately  as: 


k0~k~1Q- 


i  <  dep 

2  Cre-ar(z-z)  —  dz 


l 


dz' 


m  i,  k 

k-rA^v\  r  z 


Y  W)r  - 

L  \Az)  Jk-i 


_ar(mz_z,)dz, 


hi  ‘W 


‘z 


aT 


J  k- 1„ 


(A23) 


=  2I- 

k=  1 


Cr  k(Aep) 


,-ar(mz-fcz)  _  -ar(mz-k"1z) 


=  5^(mz) 

Equations  A22  and  A23  can  now  be  combined  to  give 


n 

s = y  ms« 

r= 1 


(A24) 
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the  expression  for  ms(r>  given  by  Equation  A23  can  be  further  simplified  as 


771  j. 

ms(r)  =  2  V  °r  (Af)  r  -ar(mz -  *z)  _  -ar(™  z-^zM 

k(Az)  i 

= 2  y Cr  k(JgP)  \e~ar(m~iz -  y  -  e-^(m_iz-fc'iz) 

A,ar  k(Az )  L 

|  2  4  k(Aen 

ar  k(Az ) 

=  m-1sMe-arm(Az)  +  2  —  _  e-arm( Az)) 

(Az) 


e-arm(Az) 


(A25) 


In  Equation  A25,  for  the  case  of  ar  =  0,  L’Hospital’s  rule  can  be  applied  to  calculate  ms^v> . 
Looking  at  Equations  A24  and  A25,  it  can  be  concluded  that  once  771  (A z)  and  777 (A  ep)  are 
known,  can  be  calculated  by  using  m~1s^r\  the  value  at  the  end  of  the  previous  loading 

step.  Thus  only  the  value  of  s(-r^  at  the  end  of  the  previous  step  needs  to  be  stored,  which  greatly 
simplifies  the  numerical  procedure. 


From  the  point  of  view  of  computer  simulation,  the  computer  code  for  the  model  should  be  able 
to  take  the  increment  in  value  of  any  six  components  of  the  stress  and  strain  tensors  as  input  and 
provide  as  output  the  incremental  values  of  the  remaining  six  components,  plastic  strain  tensor 
components,  intrinsic  time  scale  and  measure,  and  deviatoric  stress  components.  In  order  to  write 
such  a  code,  incremental  equations  are  needed.  First,  the  incremental  equation  for  the  deviatoric 
stress  can  be  given  by 


r= 1 


l(As)  =  ms  -  m'1s  =  ^  ms(r)  -  ^  m-1s(r) 

r_i  r=l 

n  n 

—  y  m-ls(r)  e-arm(Az)  +  2  _ - _ y  —  ( \  —  e-ar7”(Az)'\ 

Z_i  m(Az)  Z_i  ar  y  J 

r= 1  r—1 

n 

_  ^  m-ls(r) 

n 

=>  m(As)  =  2  "W) 


\AeP) 


l(Az)  * _ . 

r=l 
n 


_|_  m-ls(r)  ^e-arm(Az)  _  -Q 


r= 1 


where 


(A26) 


7(As)  =  771  (Act)  -  -tr(7n(A<r))/ 
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(A27) 


The  incremental  equations  for  total  stress  and  total  strain  are  given  by: 


m(Ae)  -  m(Aee)  +  m(AeP)  (A28) 

and 

m(A<r)  =  A  tr(m(Aee))  +  2/r(m(Aee))  (A29) 

The  increments  in  intrinsic  time  scale  771 (Az)  and  intrinsic  time  measure  m(A()  are  related  as: 


l(Az)  = 


l(A0 


/  (m- + rc ao,  %) 

where  0  <  /?  <  0.5  for  hardening  materials  (will  be  explained  later),  and 


(A30) 


771  (A^)  =  ^m(AeP):m(AeP)  (A31) 

Next,  a  numerical  procedure  will  be  formulated  using  the  above  incremental  equations.  Using 
Equations  A28  and  A29,  we  have 

m{Aa}  =  C[m{Ae]  -  m{A eP}]  (A32) 

where  m{Acr}is  a  vector  of  six  independent  stress  components 

m{A<j}  =  [A(T^^A(T22  A(J33  A(T23  A(Ti3  A(T^2]^  (A33) 

Vectors  m(Ae]  and  m{Aep}  are  defined  in  the  same  manner  for  total  strain  and  plastic  strain 
respectively.  C  is  a  matrix  of  elastic  constants  given  by 

'A  +  2/r  A  A  0  0 

A  A  +  2\i  A  0  0 

_ _  A  A  A  +  2\i  0  0 

L  =  0  0  0  ^0 

0  0  0  0  /r 

_  0  0  0  0  0 


0 

0 

0 

0 

0 

fU 


(A34) 


Then,  using  Equations  A27,  A28,  and  A29,  and  the  condition  of  plastic  incompressibility,  we  get 


m 


(Ac) 


2(i 


7(As)  +  m(A  eP) 
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(A35) 


where  m(Ae)  is  the  matrix  of  deviatoric  strain  increments  given  by 


m(Ae)  =  m(Ae)  -  -  tr(m(A e))I 


(A36) 


Substituting  the  expression  for  m(As)  from  Equation  A35  into  the  incremental  equation  for 
m(As),  Equation  A26,  we  have 


2p[m(Ae)  -  m(Aep)] 


=  2 


n 

1  v'1  C 


l(Az)  , _ , 

r= 1 


l(Aep) 


(A37) 


+  m_15W  ^e-arm(Az)  _ 


r= 1 


Rearranging  the  above  equation,  we  have 


M(m( Ae))  +  iy  m5(r)  (e  — arm(Az)  _ 
r=l 


’(Az))' 


r=l 


l(Az) 


Taking  the  inner  product  of  Equation  A38  with  itself,  we  get 

m(AQ2 


(A,A)  =  B‘ 


m(Az)2 


where 


A  =  p(m(Ae ))  +  ^  Y  ms(r)  (e~“rm(Az)  -  1) 

r= 1 


and 


l(Aep)  (A38) 


(A39) 


(A40) 


B=  turz)+^(l-e-«rmM) 


r= 1 


The  incremental  relation  between  z  and  (  can  be  expressed  as 


A  z  = 


(A41) 


(A42) 
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where  /?  >  1/2.  Using  Equation  A42  in  Equation  A39,  we  get 


R(Az)  =  (A,  A)  -  B2f2  +  f]A =  0  (A43) 

Equation  43,  being  a  nonlinear  equation  in  single  variable,  can  now  be  solved  for  Az  by  a  fail¬ 
safe  method  (a  combination  of  bisection  and  Newton-Raphson’s  method).  The  expressions  for 
the  derivatives  needed  to  solve  the  equation  are 


dR 

dAz 


r  I  dA 


dAz 


\  ( dB  df  \ 


(A44) 


9A  =-Y  arms(r)e"^Az 


dAz  2  Z_i  r 

r= 1 


(A45) 


dB 

dZ( 


ii 

=  li- 1-  ^  Cr  e~ 


r= 1 


df(X  +  W,§)_  + 


dAz 


(A46) 


(A47) 


Once  Az  is  obtained,  it  can  be  used  to  calculate  the  desired  plastic  strain  increment 


e 


p 

ij 


AijAz 

B 


(A48) 


The  algorithm  is  depicted  in  Figure  A6.  The  main  aspect  of  this  numerical  scheme  is  that  a 
positive  root  is  guaranteed  at  every  loading  step  regardless  of  the  step  size. 
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Assign  desired  increment  to 
<r,i  or  ca 


Calculate  the  unknowns 
by  using  eq.  (IS) 


Calculate  m(Az)  by  applying 
fail-safe  method  to  cq.  (29), 
using  eqs.  (30-33)  to  calculate 
the  slope.  Then  calculate  the 
new  value  of  '“( AcJV  )by 
eq.  (34) 


Calculate*" (A:)  and  ”‘ ( A£ )  using  eqs. 
(16)  and  (17)  resp.,  and  the  remaining 
unknowns  by  eq.  (19) 


Jt)  T  (AlTy) 

+  '"I 


•£.,  -m  ,c.,+m(AeiJ) 


m  ,  _  m  —  1 


z  +  m(Ar) 


Also  update  by  using  eq.  (12) 


Figure  A6.  Flowchart  of  the  algorithm  used  for  the  endochronic  numerical  implementation 
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A5.2.  Calibration  of  the  Material  Kernel  and  Hardening  Functions 
A.5.2.1.  Calibration  of  the  Kernel  Function 

As  discussed  previously,  the  endochronic  kernel  function  p(z)  is  defined  in  the  form  of  a  sum  of 
a  finite  number  of  exponentially  decaying  functions  of  the  form: 


f  L 

p(o = y 


Cr  e" 


r=l 


(A49) 


In  order  to  calibrate  this  function,  we  need  to  find  the  value  of  parameters  Cr  and  ar  in  such  a 
way  that  the  difference  between  the  model  predictions  and  experimental  measurements  is 
minimized.  For  most  practical  purposes,  the  energetic  composite  materials  should  be  subjected  to 
axial  loading,  hence  calibration  with  respect  to  uni-axial  cyclic  loading  experimental  data  should 
be  sufficient.  According  to  Jao  et  al.  [A10],  an  optimization  problem  can  be  formulated  in  the 
following  manner 


m 


subject  to 

min  F  (v)  —  ^  wL  [er( v,  eexp)  -  aexpf 

i= 1 

(A50) 

>4^ 

II 

p 

ci. 

ii 

p 

K> 

"S 

(A51) 

gi(v )  <  0 ;Z  =  1,2... q 

(A52) 

where  F(v)  is  the  objective  function  to  be  minimized,  uy  represents  the  weight  to  be  applied  to 
ith  data  point,  v  =  {Cr,ar}T  is  a  vector  of  the  unknown  kernel  parameters,  a(v,s)  and  <rexp 
represent  the  model  prediction  and  experimental  measurement  of  stress  respectively,  and  fj  and 
gt  respectively  represent  the  equality  and  inequality  constraints  enforced  on  the  kernel 
parameters.  Unlike  Jao  [A10],  it  has  been  assumed  here  that  the  hardening  function  is  already 
calibrated  and  known,  and  hence  the  only  parameters  that  need  to  be  determined  are  the  kernel 
parameters. 

The  above  minimization  problem  can  be  solved  efficiently  in  MATLAB  by  using  the  inbuilt 
constrained  minimization  function  fmincon.  The  complete  procedure  can  be  summarized  in  the 
following  steps: 

•  Choose  the  number  of  kernel  terms  and  take  an  initial  guess  v0  of  kernel  parameters. 

•  Specify  the  upper  and  lower  bounds  on  the  kernel  parameters.  These  are  the  constraints 
gi(v)  imposed  on  the  objective  function  (There  are  no  equality  constraints  in  this 
particular  problem,  hence  j  =  0). 

•  For  each  loading  step,  calculate  the  uni-axial  stress  from  the  3D  endochronic  numerical 
implementation,  using  uni-axial  strain  input  from  experimental  data  to  obtain  er( va,  £exp). 

•  Obtain  the  kernel  parameters  that  minimize  the  objective  function  F(v)  using  fmincon. 
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A5.2.2.  Calibration  of  the  Hardening  Function 

After  determination  of  kernel  parameters  Cr  and  ar,  the  hardening  function  is  calibrated 

by  using  the  transient  cyclic  hardening  response.  A  large  number  of  data  points  {.(J11,£11)  are 
read  from  the  experimental  cyclic  hardening  curves  at  different  strain  rates.  Then,  stress  and 
strain  increments  {Aollt  A£X1)  are  determined  between  two  consecutive  points  for  the  entire  data. 
Ae^l  is  then  determined  for  each  of  the  values  by  using  Equation  A32.  Now,  by  substituting 
771  (A z)  in  Equation  A39  using  Equation  A30,  we  have 


»nm(A 

=  3 
3 

+  ~ 


x(a4) 


(A53) 


where  m(AQ  —  ^  ^(Ae^)).  Since  171  (Ae^)  and  7n(Acr)11  are  already  known  for  all  the  data 

points  and  m-1s ^  is  determined  in  the  previous  step,  Equation  (44)  reduces  to  a  nonlinear 
equation  in  only  one  unknown,  /.  Thus  the  value  of  /  can  be  determined  at  each  step  for 
complete  experimental  data  at  each  strain  rate  by  numerically  solving  Equation  (44)  by  either  the 
bisection  method  or  Newton-Raphson’s  method.  It  can  be  assumed  that  for  a  particular  step 
(m  —  1 )  ->  m,  the  calculated  value  of  /  corresponds  to  a  value  of 

m(  =  +  /?(m(AO)  (A54) 

where  /?  =  1/2  is  a  reasonable  approximation  (Hsu  et  al.  [A9]).  Now,  values  of  /  corresponding 
to  various  data  points  on  the  stress-strain  curves  of  each  strain  rate  can  be  plotted  against  (£,  <f), 
and  the  hardening  function  /(£,£)  is  determined  by  2D  linear  piecewise  fitting  to  the  plot. 

A.5.2.  Staggered  Algorithm  for  Simultaneous  Calibration  of  the  Kernel  and  Hardening 
Functions 

In  the  previous  sections  we  have  discussed  the  calibration  procedure  for  the  kernel  and  hardening 
function  parameters,  and  have  seen  that  the  two  functions  are  interdependent  for  the  purpose  of 
calibration.  We  assumed  there  that  the  calibrated  hardening  function  was  available  for  the  kernel 
function  and  the  calibrated  kernel  function  was  available  for  the  hardening  function.  However,  in 
practice  we  would  need  to  calibrate  both  of  the  functions  at  the  same  time.  This  problem  was 
efficiently  tackled  through  a  staggered  algorithm  that  could  simultaneously  calibrate  both  the 
kernel  and  hardening  functions.  Figure  A7  provides  a  description  of  the  algorithm.  The  iterative 
procedure  starts  with  an  initial  guess  of  the  kernel  parameters,  which  along  with  uni-axial  cyclic 
loading  experimental  data,  is  provided  as  an  input  for  calibration  of  the  hardening  function.  The 
calibrated  hardening  function,  along  with  the  initial  guess  of  kernel  parameters  and  experimental 
data  serve  as  an  input  to  the  optimization  problem  that  yields  a  new  set  of  kernel  parameters.  The 
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procedure  is  iterated  until  the  error,  denoted  by  the  Z2  norm  of  the  difference  between  the  kernel 
parameters  of  the  current  and  previous  guess  converges  within  a  specified  tolerance. 


Figure  A7.  The  staggered  algorithm  used  for  the  simultaneous  calibration  of  the 
endochronic  kernel  and  hardening  functions 
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